#!/bin/sh

BLAST_HOME=/Users/david/Tools/ncbi-blast-2.2.25+
MCL_HOME=/Users/david/Tools/mcl
PATH=$PATH:$BLAST_HOME/bin:$MCL_HOME/bin

sample=$1
reference=$2
pre=$3
inflation=$4

echo "####################################################" > process.log
echo "#         Date/Time : \c" >> process.log
date >> process.log
echo "#             Sample: $sample" >> process.log
echo "#          Reference: $reference" >> process.log
echo "#  MCL Pre-Inflation: $pre" >> process.log
echo "#      MCL Inflation: $inflation" >> process.log
echo "####################################################" >> process.log

sample_name=`echo $sample| cut -d'.' -f 1`

thisTime=0
startTime=0
stopTime=0
elapsed=0

setTime() {
  set `perl -MTime::HiRes=gettimeofday -e 'print int (1000*gettimeofday()).qq(\n);'`
  thisTime=$1
}

startTime() {
  setTime
  startTime=$thisTime
}

stopTime() {
  setTime
  stopTime=$thisTime
  elapsed="$(expr $stopTime - $startTime)"
  echo $1: $elapsed ms | tee -a timings.log
  startTime
}

rm -rf out > /dev/null
mkdir out
rm -rf db > /dev/null
mkdir db 

startTime

# [A1]
cat $sample | awk '{print ">" substr($0,2);getline;print;getline;getline}' > out/$sample_name.fa
stopTime "[A1] FASTQ => FASTA Conversion"

makeblastdb -dbtype nucl -in out/$sample_name.fa -out db/$sample_name >> process.log
stopTime "[A1] Make Sample DB"

# [A2]
blastn -word_size 11 -reward 2 -penalty -3 -gapopen 5 -gapextend 2 \
  -query $reference \
  -db db/$sample_name \
  -outfmt '6 sseqid qseq' \
  | awk '{print ">" $1;print $2;getline}' > out/intermediate.fa
stopTime "[A2] Query"

# [A3]
makeblastdb -dbtype nucl -in out/intermediate.fa -out db/intermediate >> process.log
stopTime "[A3] Make Intermediate DB"

# [A4] 
rm -rf tmp > /dev/null
mkdir tmp
cd tmp
split -a 7 -l 2 ../out/intermediate.fa
cd ..

rm -rf data > /dev/null
mkdir data
for f in `ls tmp/*`
do
  set `cat $f | grep '>' | awk '{print substr($0,2)}'`
  name=$1
  blastn -word_size 11 -reward 2 -penalty -3 -gapopen 5 -gapextend 2 \
    -query $f \
    -db db/intermediate \
    -outfmt '6 sseqid evalue' \
    -out data/$name
done
stopTime "[A4] Query"

rm -rf output > /dev/null
mkdir output
java -cp ../../../target/classes com.attrisoft.sanger.wrapper.Main
stopTime "[A5] Matrix Conversion"

# [A6]
cd out
mcxload -abc data.abc --stream-mirror -write-tab data.tab -o data.mci >> ../process.log 2>&1
mcl data.mci -I $inflation -pi $pre -o out.data.mci >> ../process.log 2>&1
mcxdump -icl out.data.mci -tabr data.tab -o dump.data.mci >> ../process.log 2>&1
cd ..
stopTime "[A6] MCL Execution"

java -cp ../../../target/classes com.attrisoft.sanger.wrapper.Clusters
stopTime "[A6] Cluster Extraction"

cp out/intermediate.fa output
mv timings.log output
mv process.log output

rm -rf input > /dev/null
mkdir input
mv $sample input
mv $reference input
